arXiv: 1506.030lOvl [astro-ph.GA] 9Jun2015 


Draft version June 10, 2015 

Preprint typeset using DTfTjX style emulateapj v. 5/2/11 


GRAVITATIONAL ENCOUNTERS AND THE EVOLUTION OF GALACTIC NUCLEI. II. 

CLASSICAL AND RESONANT RELAXATION 

David Merritt 

Department of Physics and Center for Computational Relativity and Gravitation, Rochester Institute of Technology, Rochester, NY 

14623 

Draft version June 10, 2015 

ABSTRACT 

Direct numerical integrations of the Fokker-Planck equation in energy-angular momentum space are 
carried out for stars orbiting a supermassive black hole (SBH) at the center of a galaxy. The algorithm, 
which was described in detail in an earlier paper, includes diffusion coefficients that describe the 
effects of both random (“classical”) and correlated (“resonant”) encounters. Steady-state solutions 
are similar to the Bahcall-Wolf solution, n(r) oc r -7 / 4 , but are modified at small radii due to the 
higher rate of diffusion in angular momentum, which results in a low-density core. The core radius 
is a few percent of the influence radius of the SBH. The corresponding phase-space density f{E,L) 
drops nearly to zero at low energies, implying almost no stars on tightly-bound orbits about the 
SBH. Steady-state rates of stellar disruption are presented, and a simple analytic expression is found 
that reproduces the numerical feeding rates with good accuracy. The distribution of periapsides of 
disrupted stars is also computed. Time-dependent solutions, f(E,L,t), are also computed, starting 
from initial conditions similar to those produced by a binary SBH. In these models, feeding rates 
evolve on two timescales: rapid evolution during which the region evacuated by the massive binary is 
refilled by angular-momentum diffusion; and slower evolution as diffusion in energy causes the density 
profile at large radii to attain the Bahcall-Wolf form. 


1. INTRODUCTION 

Paper I in this series (jMerrittJ 120151 ) presented a numerical algorithm for integrating the Fokker-Planck equation 
describing f(E, L, t), the phase-space density of stars orbit ing a S B H at the center of a galaxy. The algorithm described 
in Paper I was similar to that in the pioneering study of iCohn fc Kulsrud ( 197 8). but with a few modifications. Loss 
of stars into the SBH was treated more carefully, by adopting a logarithmic grid in angular momentum and by 
incorporating a more precise expression for the loss-cone flux. In addition, the diffusion coefficients in energy, E , 
and angular momentum, L , were allowed to have more general forms than those derived in the classical theory of 
Chandrasekhar, Henon, Spitz er an d others, all of whom assumed random (uncorrelated) interactions, and (with a few 
notable exceptions, e.g. iLeel ( 1969 11 ignored relativistic corrections to the equations of motion. 

Two characteristic length scales are commonly associated with a supermassive black hole (SBH) at the center of a 
galaxy. The gravitational radius r g , 
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is the length scale set by Einstein’s equations for a relativistically compact object. The (gravitational) influence radius, 
rinfi, has two standard definitions: either r; n fl = 77 ,, where 
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or Hnfl = r mi defined implicitly via 


M*(r < r m ) = 2 M.. 


( 3 ) 


The first of these is expressed in terms of a, the one-dimensional velocity dispersion of stars in the galactic nucleus, 
while the second is defined as the radius containing a stellar (distributed) mass equal to twice M.. In the Milky Way, 
fh ~ Tm ~ 10°' 5 pc (|Schodel et al.ll2009i : iChatzopoulos et alJl2015h . 

The classical theory of gravitational encounters is valid at distances r > r m from a SBH. In this regime, random 
encounters imply a similar timescale for changes in both E and L: the “two-body” relaxation time T r , given by 


T r = 


0.34 cr 3 
G 2 m+p In A 


1.2 x 10 


10 


100 kms 


-1 


-1 


10 5 M Q pc- 3 7 \m q J V 15 


m* \ 1 / In A 


-1 


yr- 


( 4 ) 


Here p is the stellar mass density, m* is the mass of a single star, and In A is the Coulomb logarithm dChandrasekliad 
11942H . In Figure 1, the region where equation d4| defines the timescale associated with gravitational encounters is 
labelled “Newton.” 
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If one imagines approaching ever more closely to the SBH, the unperturbed orbits change in well-defined ways, 
implying corresponding changes in the dominant mode of gravitational interaction between stars. Starting at a radius 
of - io-V^B orbits are so nearly Kepl erian that th e ass umption of uncorrelated encounters breaks down. This is the 
regime of “resonant relaxation” (iRauclr fo Tremainej|l99'6 ) in which changes in L can occur on much shorter timescales 
than changes in E. The corresponding spatial region is labelled “Kepler” in Figure 1. 

Still closer to the SBH, the lowest-order effects of general relativity (GR) begin to make themselves felt. Orbits 
experience planar (apsidal) precession due to the 1st post-Newtonian (1PN) corrections to the equations of motion. 
One consequence is that the coherence time for the resonant interactions described above is determined by GR in this 
region. In addition, the most eccentric orbits at a given energy will precess due to GR at a higher rate than most 
other orbits of simil ar ener g y, causing the former to behave in qualitatively different ways than the latter in response 
to perturbations (iMerritt et alJ f201~Tl l. The region where gravitational encounters are strongly affected by these 1PN 
relativistic effects is labelled “Schwarzschild” in Figure 1; this region has an outer radius of — 10 -2 r m . 



Fig. 1 - Sketch of the different regions 
around a massive black hole at the center 
of a spherical galaxy. Each region is de¬ 
fined in terms of the dominant mechanism 
by which gravitational encounters change 
the orbits of stars. The outer circle is the 
black hole’s gravitational influence radius 
r infl an d the inner circle is its gravitational 
radius r g . 


Another change in the character of the unperturbed motion occurs still closer to the SBH, where relativistic frame- 
dragging causes orbits to precess nodally about the spin axis of the SBH. While the importance of frame dragging 
for gravitational encounters has hardly begun to be explored, one consequence is known: within a certain distance, 
Lense-Thirring torques from a Kerr black hole dominate the Ne wtonian torques that would otherwise (via “ vector 
resonant relaxation”) be responsible for changes in orbital planes (IMerritt et, al.l[2?Tl0l : llVlerritt fc Vasilievl 1201^1 . This 
region, of outer radius — 10 -3 r m , is labelled “Kerr” on Figure 1. 

At still smaller radii, PN terms higher than 2nd order can become important, implying changes in E and L due 
to gravitational-wave emission. Within a distance of perhaps 10 2 r s from the SBH, relativistic corrections can not be 
adequately treated via a Newtonian or post-Newtonian formalism. This regime is labelled “Einstein” in Figure 1. Of 
course the number of stars (or compact objects) present at any time in this region is likely to be small. 

When evaluating the response of stellar orbits near a SBH to gravitational encounters, most researchers have applied 
the classical expressions for the diffusio n coefficients , implicitly assuming that those expressions remain valid arbitrarily 
close to the SBH. An example is the lBahcall & Wolf (1976!) steady-state solution for a single-component stellar cluster: 

/ oc |£| 1/4 , n oc r" 7/4 (5) 

which was derived using the iHenonl ( 1961 .) expressions for the orbit-averaged diffusion coefficients (A E), ((A E) 2 ). 
This approach is defensible given that the form taken by the diffusion coefficients in the regions below “Newton” in 
Figure 1 can not yet be derived fro m first p ri nciples. __ __ 

However, recent numerical work (|Hamers. Porte gics Zw art fe Merrittl 12014 : Merritt 2015|) has yielded approximate 
and fairly general expressions for the angular momentum diffusion coefficients in the “Kepler” regime. Those ex- 

1 More precise estimates of radii like this one are presented below. 
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pressions can be included in a Fokker-Planck description and used to evolve f{E,L ), yielding solutions that are 
valid—if not at all ra dii— at least to distances ~te n times closer to the SBH than classical solutions like those of 
iBahcall fe Woli (I1976T) and I Cohn fe Kulsrudl (119781 ). The diffusion coefficients adopted in the present paper are in 
fact valid even into the “Schwarzschild” regime of Figure 1, in the sense that they correctly account for the effects of 
1PN apsidal precession on the coherence time. However, no attempt is made here to treat the effects of “ano ma lous 
relaxation,” the qualitatively different way in which low-T orbits evolve in the Schwarzschild regime (jMerritt et al.1 
120111: lHamers. Porteeies Zwart fe Merrittll2014D . 

Section [2] reviews the numerical algorithm used here; further details are given in Paper I. In § [3] timescales associated 
with gravitational encounters in the “Newton,” “Kepler” and “Schwarzschild” regimes are derived for the case of stars 
in a Bahcall-Wolf cusp around a SBH. Section [4] presents steady-state and time-dependent solutions for f{E,L). 
Section [5] discusses some implications of the results obtained here for real stellar systems and jjGjsums up. 


2. EVOLUTION EQUATIONS 

The numerical algorithm for integrating the orbit-averaged Fokker-Planck equation is described in detail in Paper 
I. Features of the algorithm that are most relevant to the current study are reviewed here. 

Stars are assumed to have a single mass, ?n*, and to be close enough to the black hole (SBH) that the gravitational 
potential defining their unperturbed orbits is 

<F(r) = — = —i/j{r) (6) 

r 

with M m the SBH mass, assumed constant in time. Unperturbed orbits respect t he two is olating integrals E , the 
energy per unit mass, and L, the angular momentum per unit mass. Following iCohn fc Kulsrudl ( 1978 ) these are 
replaced by £ and TZ where 


v 2 L 2 

£=-E=- Y +^(r), 


(7) 


L c (£) is the angular momentum of a circular orbit of energy £ so that 0 < TZ < 1. £ and 7 Z are related to the 
semimajor axis a and eccentricity e of the Kepler orbit via 


The orbital (Kepler) period is 


GM . 2 t -p 
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2S 
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and L c = GM m /y/2£ = y/GM m a. Spin of the SBH is ignored. 

The time dependence of the phase-space number density of stars, f{£, 1Z), is described by the orbit-averaged Fokker- 
Planck equation 


d f d d 

d f df d f df 

~<t>£ = D££— + De-ji—— + Dgf, —<t>Ti = Dug-— + D'lm—— + D-jif ( 10 ) 

o£ oTZ a£ olZ 

with flux coefficients 

De = -m - j(«A £f) + —mf) + ~(A£AK) , 

D n = -(ATI) - y (A£A 1Z) + l-^(A£ATZ) + ^<(A 1Z) 2 ) , 

D ££ = \{{A£) 2 ) , D £n = D ns = \{A£ATZ) , D nn = \({A1Z) 2 ) (11) 


and J = \/2tt 3 G 3 M, 3 £ 5 / 2 (lMerrittll2013l 5.5.1). Quantities in ( ) are orbit-averaged diffusion coefficients, which are 
expressed as 


{AS) = (A£)ck, (( AS) 2 ) = ((A£) 2 ) C k, {AS ATZ) = {A£A1Z) C k, 

(A 1Z) = (A72.)ck + (A^) rr , {{ATZ) 2 ) = {{ATZ) 2 ) C k + ((A^) 2 ) RR . (12) 

The subscript CK indicates that the diffusion coefficient is computed as in Cohn fc Kuls rudl ( T978D : their derivation 
was based on standard assum ptions about randomness of encounters (Roscnb luth et al.lll957D . The subscript RR 
refers to “resonant relaxation” dRauch fc Trem aine 1993). The resonant diffusion coefficients are assumed to have the 
simple, separable forms 


(A^)rr = 2A{£) (1 - 2 TZ ), ((A^) 2 ) rr = AA{£)TZ{1 - TV). 


(13) 
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The term containing the £ dependence is 


A(a) = 
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Here IV = /V(r < a) is the number of stars instantaneously at radii smaller than a, M* = m+N, P is the Kepler 
(radial) period, and t co h is the coherence time, defined as 


f - 1 =f - 1 Mf - 1 
b coh ~ l, coh,M ^ Koh,S 


M. 

tcoh,M(a) — P 
Nm+ 


1 a 

^coh,s(&) = T77 P- 
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(15) 


f co h,M is the mean precession time for stars of semimajor axis a due to the distributed mass around the SBH (“mass 
precession”), and t C oh,s is the mean precession time due to the 1PN corrections to the Newtonian equations of motion 
(“Schwarzschild precession”). 



Fig. 2. — Angular momentum diffusion coefficients plotted as functions of 72. = L 2 /L 2 . Open (red) circles show equations 1131 . the adopted 
expressions for the resonant diffusion coefficients. Lines show the Cohn-Kulsrud diffusion coefficients, in models having n(r) oc r 1 , i.e. 
/ oc and 7 = {1/2, 1, 3/2, 2}. Curves on the left are normalized to the same y— value at 72. = 0, and curves on the right are 

normalized so as to have the same peak value. In the left panel, curves with largest 7 have the smallest ordinate value at 72 = 1; in the 
right panel, the curves with 7 = 0.5 and 7 = 2 have the largest ordinate values at 72 ~ 0.05, while the two curves with 7 = 1 and 7 = 3/2 
lie slightly below. 


The functional forms chosen for the resonant diffusion c oefficients in equations (113l) - (115l) wer e show n in Paper I to 
reproduce the numerically-extracted diffusion coefficients of lHamers. Portegies Zwart fc Merritt ( 201 1 i. at least in the 
particular set of scale-free nuclear models considered by those authors. Figure [2] makes another comparison: between 
the ^-dependence of the resonant diffusion coefficients, equations (THU) , a nd t he ^-dependence of the classical diffusion 
coefficients. The latter were computed from equations (24)-(25) of lCohn fe Kulsrudl (119781) . assuming scale-free forms 
for / and n: 

f oc £ 7-3 / 2 , n oc r -7 (16) 

and a 1/r potential. The different curves in Figure [2] were normalized as described in the figure caption; of course 
the classical and resonant diffusion coefficients can have very different amplitudes. It is remarkable that the classical 
diffusion coefficients are very well fit by the same simple functions of 1Z that were adopted for the resonant diffusion 
coefficients; and furthermore that there is so little dependence of the former on 7 . These results may lend an extra 
degree of confidence to the functional forms assumed here for the resonant diffusion coefficients. 

Loss of stars into the SBH is controlled by the choice of n c , the radius of the physical loss sphere around the SBH, 
and by the conditions imposed on / at the loss-cone boundary, 1Z = 1Z\ C (£), defined as 

k '‘< £ >= 2 IH£)’ < ir > 

7Z\c is the normalized angular momentum of an orbit with (Newtonian) periapsis at r\ c . The /^.-directed flux of stars 
across the loss-cone boundary is 


F{£) d£ = -J{£)MU lc ) d£ = -J(£) </> k , i c (£) d£. 


( 18 ) 
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Two quantities that play important roles in angular momentum diffusion near the loss-cone boundary are T >, 

ftlc) 


V{£) = ^ 

v ’ 277 


n=n lc 


72-ic 


(19) 


and , 


qi c(£) = 


P{£) V{£) 

Kte{£) 


( 20 ) 


T>~ 1 is effectively an orbit-averaged, angular momentum relaxation time at energy £. The quantity qi c measures the 
change in angular momentum per orbital period, compared with the size of the loss cone. The loss-cone boundary 
conditions adopted in all the integrations presented here were the “Cohn-Kulsrud boundary conditions” defined in 
Paper I. No attempt is made to solve for / inside the loss cone, i.e. at 77 < 77i c , since / does not satisfy Jeans’s 
theorem in this region. 

Solutions are obtained on a (N x x N z ) grid in (X, Z), where 


X = In R = In 



Z = ln(l+/3£*) = In (1 +/3£/c 2 ) . 


Integrations presented here used N x = N z 
The code adopts units such that 


64 grid points. 

G = M,=c= 1 


( 21 ) 

( 22 ) 


allowing the results to be scaled to different masses of the SBH. Dimensionless parameters that must be specified 
before the start of an integration include m*/M,, In A and 0i c = r\ c /r g . 

It is important to emphasize that all the results presented in this paper assume a single mass for the stars. One 
reason for this simplification is the current, poor state of knowledge about the form of the resonant diffusion coefficients 
in systems containing a range of stellar masses. The mass dependence of the classical diffusion coefficients is of course 
known; it implies a rate of energy loss for massive objects that scales in proportion to their mass, leading to segregation 
of the more massive objects tow ard the galaxy center. In steady-state models of t he Milky Way nucleus that contain 
a realistic stellar mass function (Frcitag ct al. 20061: iHonman fc Alexand er 2006a]), mass segregation implies that the 
total density is dominated by the heaviest stellar remnants, ~ 10 Mq black holes, inside a sphere of radius ~ a few 
xl 0~ 3 pc. 


3. IMPORTANT QUANTITIES IN A BAHCALL-WOLF CUSP 

The steady-state numerical solutions presented later in this paper can be described as modifications of the classical 
Balicall-Wolf solution, equation ©. Here we evaluate some important quantities associated with angular momentum 
diffusion in a nucleus with 

p(r) = m*n(r) oc r _7//4 , f(£) oc S 1 ^ 4 , i>{r) = ——(23) 

r 

a unmodified Bahcall-Wolf (1976) cusp. Like Bahcall and Wolf, we ignore here the L— dependence that a realistic / 
would necessarily have due to capture by the hole. 

Define r m as the radius containing a mass in stars of 2M,. The number density is 


n(r) = n 0 



5 M. 1 f r \~ 7/4 
87 r m* r^ \r m ) 


(24) 


The number of stars with instantaneous radii less than r, or semimajor axes less than a, are given respectively by 
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The phase-space number density is 
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and the distribution of energies is 


V2ti 


N(£) d£ = 4ir 2 p{£)f(£) d£, p{£) = ( GM.) 3 £~ 5 / 2 . 


The mass coherence time, equation (1151) . is 
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and the Schwarzschild coherence time, equation (115b '). is 
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Figure [ 3 ] plots these times, as well as the coherence time f co ^ = t co ^ M + t c ° h s defined in equation m- 


^coh(^) — \ 
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(28) 

(29) 
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(31) 

(32) 


assuming M. = 4 x 10 6 M ( ? tl and for two choices of r m : 1 pc an d 10 pc, which probably bracket the actual value at the 
Galactic center (ISchodel et al.l[20091 : IChatzopoulos et al.ll2015h . 




Fig. 3. — Characteristic times inancc ,/4 nucleus. Parameters are M. = 4 x 10 6 Mq, m* = l.OMg, r m = 1 pc (left) and r m = 10 pc 
(right); the value of r m at the Galactic center probably lies between these two values. The curve labelled ri c assumes a capture radius of 8r g , 
appropriate for compact objects; tidal disruption of Solar-mass stars would occur at gr eate r dist ance s. The two solid (black) curves labelled 
t co h show the coherence times due to mass- and Schwarzschild precession, equations (1281 and unii respectively; the mass (Schwarzschild) 
coherence time has the smaller (larger) value at large radii. Dot-dashed (red) line is 10 5 times the Kepler period P. Thick (magenta) 
curve sho ws t he overall coherence time, equation 102}. The three curves labelled T> 1 are equations rm . i36ii . and 03: magenta curve is 
equation 09}. In the case of the curves showing timescales, “radius” means “semimajor axis.” 


The quantity T> defined in equation (1191) is effectively an inverse, orbit-averaged, angular momentum relaxation time. 
In the case of classical relaxation, the CK diffusion coefficients imply, in the limit 1Z —> 0, 
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where T = 47r(Gm*) 2 In A, i.e. 


^ 5C r(ll/4) G 2 m*M. In A c<1/4 

/^— TV C/A\ _,7/4 S/4 C 1 


r(5/4) ( GM ,) 7/4 ri /4 


2? 


-1 , 


: 1.61 X 10 9 


M. 


1/2 


5/4 


M q 


/In A 


fe) 


1/4 


4 x 10 6 A2© y y 3 pc 

In the case that diffusion is dominated by resonant relaxation, equations m and ED imply 
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Setting t coh = t C oh,M gives 
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Equating (1341) and (1361) (which assumes t co h = t C oh,M) yields an estimate of the radius below which resonant relaxation 
dominates classical relaxation: 
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This radius can be identified with the sphere labelled “Kepler” in Figure 1. 

Adopting equation (1151) for the overall coherence time, we can write an expression that is valid throughout the 
resonant-relaxation-dominated regime: 
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The diffusion time associated with resonant relaxation reaches a minimum when 
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(40) 


slightly smaller than the radius at which f C oh,M = icoh.s- Either of these radii can be associated with the sphere labelled 
“Schwarzschild” in Figure 1. 


4. RESULTS 

4.1. Steady-state solutions 

ICohn fe Kulsrudl (|1978T ) obtained various steady-state solutions for /(£, TZ), assuming classical relaxation, and with 
parameters chosen to represent stars orbiting a massive black hole in a globular cluster. As they noted, an algorithm 
that ignores the contribution of the distributed mass to the gravitational potential can not be expected to correctly 
represent the solution for / at low binding energies; that is, beyond the black hole’s gravitational influence radius. In 
all of their integrations, the outer boundary condition was taken to be 

f(E = 0,K) = f 0 


(41) 































and they identified /o with n/(2n(v 2 )) 3 ' 2 ; n and ( v 2 ) are respectively the number density and mean-square stellar 
velocity in the cluster core, where the density was assumed to be constant with radius. Cohn & Kulsrud interpreted 
their outer boundary condition as describing a fixed, Maxwellian velocity distribution at large distances from the black 
hole. 

Solutions in this section were computed using a similar outer boundary condition: 

f*(£ min ,K,t) = f*(£ min ,K,0). (42) 

Here, /* is the dimensionless phase-space density, and £ m in is the minimum value of £ on the energy grid0 An 
approximately equivalent statement is that the outer boundary condition consisted of specifying a fixed mass density 
at the outermost grid radius. 

Near the loss-cone boundary 1Z = TZ\ C {£) 1 the Cohn-Kulsrud conditions were imposed, in the manner described in 
detail in Paper I. 

The initial conditions for f(£,lZ) were based on an isotropic power-law model, n oc r _T , / oc <? 7-3 / 2 , but with a 


simple modification to account for the presence of the loss 

f(£,n) = o, 



cone, namely 

n<n lc {£). (43) 



Fig. 4.— Left panel: Steady-state density profiles for integrations with three values of the outer density normalization and six values of 
m*/M.. From bottom to top in each set, m* = (0.3,1, 3,10, 30, 100}Mq, assuming M, = 4.0 X 10 6 Mg. Right panel: The three curves 
from the left panel with m* = lMg have been replotted, and rescaled vertically to give the same mass density at r = 1 pc. Dotted line 
has the Bahcall-Wolf slope, p oc r~ 7 / 4 . Vertical tick marks indicate where t co h M = /: (;0 | l g. 

Some examples showing the time-evolution of n(r, t), starting from initial conditions similar to these, were presented 
in Paper I. Here we focus on the steady states. To minimize integration times, the value of 7 defining the initial 
conditions was set to 7/4, close to the expected, steady-state value. Integrations differed in their choice of two 

parameters: the outer mass density; and m*/M # . Assuming M, = 4 x lO 6 Af 0 (the code sets the SBH mass to one), 

the adopted values of to* were 

to* = {0.3,1,3,10,30,100}M Q . (44) 

In principle, one could identify each value of to* with stars of a certain type and estimate the corresponding tidal 
disruption radius. Instead, the radius ?'i c of the loss sphere was chosen to be a fixed multiple of r g in in all integrations, 
01c = n c /r g = 15, i.e. 

nc = 16r,»2.9xlO-«( Jlr ^)pc, (45) 

roughly the value of the tidal-disruption radius for a Solar-type star at the Galactic center. Many properties of the 

steady-state solutions, including the rate of loss of stars to the SBH, are expected to depend only logarithmically on 

He- 


It is likely that Cohn Sz Kulsrud also enforced their boundary condition at a finite £ m inj and not at £ = 0. 
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Figure [I] shows steady-state density profiles for integrations with three different outer boundary conditions, corre¬ 
sponding to mass densities at one parsec of roughly 

{1.9 x 10 4 ,3.5 x 10 5 ,6.1 x lO 6 }M 0 pc- 3 . (46) 

These values probably bracket the actual value in the Milky Way (|Schodel et alJl2009t IChatzopoulos et al.l lT015h The 
left panel shows solutions for each of the 18 models, i.e., six values of m* for each choice of outer density. To a good 
approximation, the form of p(r) is determined by the (mass) density normalization, independent of m*. The reason 
can be seen by comparing equations (1341) . (1361) and (1371) . which show that the rate of angular momentum diffusion 
scales with m* in the same way for both classical and resonant relaxation, if a fixed value of r m —that is, a fixed 
mass density—is assumed. The (weak) dependence of the steady-state density profile on m* is due to the fact that 
qi c , defined in equation (BUI) , is also proportional to to*. Solutions with the smallest to* approach most closely to the 
“empty-loss-cone” form, qi c <C 1, for which 


n lc (£)<n<i 

while large values of to* imply q\ c 1 and 

f(£,lZ) « const., lZic(£) < TZ < 1, 
the “full-loss-cone” solution ( :Merritdl2013l 6.1.2). 


(47) 

(48) 



In m 



e/c 2 


Fig. 5.- Left: Equilibrium phase-space density in the integration from Figure [4] with m*/M@ = 1 and with a final density at 1 pc of 
~ 3.5 x 10 5 Mq pc -3 . Greyscale is proportional to log / and the thin (yellow) curves are contours of constant /. Thick (blue) curve is the 
loss-cone boundary and thin (red) curve is 'R.o(S). The solution grid was uniform in the plotted variables {X = InTZ, Y = ln(l + /3£/c 2 )}; 
grid centers are indicated with the dots. Right: Angular-momentum-averaged distribution functions for the three, steady-state models 
from Figure [41 with m* = Mq; the middle curve corresponds to the steady-state / plotted at left. The vertical normalization of each curve 
was chosen to give a fixed value at low energies. Dotted line shows the Bahcall-Wolf solution and vertical dashed line indicates the energy 
of a circular orbit at the assumed radius of the capture sphere around the SBH. 


For mation of a “core” is an expected consequence of resonant relaxation (Hopman & Alexander 2006b: Madigan et al. 
l20ll . In the “Kepler” regime (Figure 1), diffusion in angular momentum takes place on a short timescale compared 
with diffusion in energy. As a consequence, stars in this region are scattered into the SBH in a time short compared 
with the time for the same orbits to be repopulated by (classical) energy diffusion. An estimate of the value of a below 
which resonant relaxation dominates classical relaxation was made in equation (1381) : a ~ 3 x 10~ 2 r m . Since 

T m ~ {0.2, 2.0, 20}pc (49) 

in the models of Figure [2 the value of a at transition is predicted to be ~ {6 x 10 _3 ,6 x 10~ 2 ,6 x 10 —1 } pc. It is 
reasonable to divide these values by ~ 2 to convert from a to radii. The resulting values are quite similar to the radii 
of the cores in Figure 2) 

The right panel of Figure [I] compares the steady-state density profiles in the three integrations with m* = Mq. To 
assist in the comparison, the curves have been adjusted vertically so as to have the same density at a radius of one 
parsec. 
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Depletion of / at high binding energies should eventually result in gradients with respect to E that drive a (classical) 
flux that balances the losses due to (resonant) diffusion in L. The left panels of Figures [5] and [6] provide support for this 
statement. Plotted there are the steady-state /(£, 1Z) (Figure[5]), and streamlines of the flow in (£, TZ) space (Figure[6]), 
of the model from Figure |4] with to*/Mq = 1 and with the intermediate, large-radius density. There is a remarkably 
strong depletion of f above a certain binding energy. The right panel of Figure [5] shows angular-momentum-averaged 

/’s: 

/(£) = [ f(£,TZ) dlZ (50) 

do 

for the three steady-state models from Figure [4] with m*/Mg = 1. To a good approximation, / = 0 above a certain 
£, and so the configuration-space density at small radii has the form 

n(r) ~ r -1 / 2 , (51) 

the density of a population of stars with a single energy moving in a 1/r potential. This is approximately the central 
dependence of p on r in the profiles of Figure [4] 

The thin (red) curve in the left panels of Figures [5] and [5] is the qu antity 7 £q (£), the / = 0 intercept of the Cohn- 
Kulsrud boundary-layer solution extrapolated inside the loss cone (|Merrittll2013l equation 6.65). An “empty” loss cone 
has IZq ~ TZ\ C - At low binding energies, TZq can be seen to drop below 1Z\ C , indicating that the loss cone is becoming 
progressively fuller far from the SBH. 


CO 
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Fig. 6.- Left: Streamlines of the flux f eonation 1 1 ( 11 1 in the steady-state model of Figure[5] Right: Flux of stars into the loss cone as a 
function of energy. The three sets of curves correspond to the three models from Figure [4] Each curve has been normalized vertically to 
give a peak flux of one. Dashed curves show the contribution to <fiiz from resonant relaxation. 


The right panel of Figure [G] plots the 7?.-directed flux, fa at the loss-cone boundary as a function of energy, in the 
three steady-state models of Figure 0] Corresponding to the depletion of / at large binding energies, there is a similar 
depletion in the loss-cone flux, such that (f>Ti peaks narrowly around a certain energy. The dashed curves in this figure 
show the contribution to the flux from resonant relaxation. As expected, the resonant contribution to the flux becomes 
dominant at roughly the same energy where the depletion in / occurs. 

A quantity more directly related to the loss rate than (/n,\c{£ ) is F(£) = —J{£)(j)n,\c{£)'i equation (fT51) states that 
the integral of F(£) with respect to energy yields N. The left panel of Figure [7] plots \£F{£)\ for each of the 18 
steady-state models of Figure [I] This quantity can be interpreted as the contribution to the total loss rate from stars 
in the energy interval d£/£ 1 or equivalently, da/a. It is clear from this figure that there is a substantial contribution 
to the feeding rate from stars at large radii, hence in the classical regime, particularly in the case of small to*. i.e. an 
empty loss cone. As ?n* is increased (at fixed p), the radius of transition from full- to empty loss cones drops; the 
curves peak, roughly, at this radius. In the empty-loss-cone regime, near the SBH, loss rates (measured in stars per 
year) are nearly independent of to*. This follows from the dependence T> ~ to*, noted above, and the fact that for a 
fixed mass density, the number of stars scales as to” 1 . Far from the SBH, in the full-loss-cone regime, the expressions 
derived below imply 


\4>tz,\c\ oc to* 1 £ 11/4 , \£J4>k,\c\ oc to* 1 £ 5/a , £ ->■ 0. 


(52) 
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The total loss rate diverges in the case of an everywhere-empty loss cone, to* — > 0, but is finite for finite ?n*. 

The right panel of Figure [7] shows integrated loss rates for the same set of models. As expected, N attains a well- 
defined limit at low £, i.e. large a, in each model. These values are listed in Table [l] and plotted against to* in 
Figure [S] In three of the models, the adopted energy grid probably did not extend to low enough values of £ to yield 
accurate values for the total loss rate; these numbers have been placed in parentheses in Table [0 




Q Ag Q / r g 

Fig. 7.— Local (left) and integrated (right) loss rates, defined as number of stars per year, for the equilibrium models of Figure PH 
assuming M m = 4 X 10®Mq. In each set of curves, the value of in* increases from top to bottom. The energy £ = GM,/(2a) at which 
91c = | In 77i c | is indicated by a circle in the curves on the left. 


It is useful to have an approximate analytic expression for the loss rate. Since most of the stars lost to the SBH 
are orbiting in the “Newton” regime prior to capture (Figure [7]), it is reasonable to adopt the classical expressions for 
the angular-momentum diffusion coefficients at all energies. We simplify the derivation even more by (i) adopting for 
the density profile a Bahcall-Wolf cusp, unmodified by resonant relaxation and by loss-cone effects; and (ii) assuming 
that stars at any given energy are either in the empty-loss-cone, or the full-loss-cone, regimes at the time of capture. 
In the full-loss-cone (FLC) regime, / is assumed to be independent of 1Z , f FLC (£,lZ) = /(£). I n the empty-loss-cone 
(ELC) regime, equation (ITHl for / can be written 


m,£) 


m i) 

]n(l/ftic) 


In 



m 

ln(l/7?-ic) + TZ\ C ~ 1 


In 



Hie < U < 1 


(53) 


with /(£) defined as in equation ([50]) . The differential loss rate, F(£), is defined such that the number of stars lost, 
per unit of time, from orbits with energies in the range £ to £ + d£ is F(£ )d£. In the FLC regime, the loss rate is 
equal to the orbital draining rate: 

F flc (£) = 4t w 2 L 2 c (£)1l lc (£)f(£) = P(£)~ l n lc (£)N{£) (54) 

(lMerrittll2013l . equations (6.10,6.72)). The loss rate from stars in the ELC regime is 

ELC 4t t 2 L 2 c (£)P(£)V(£) 7 V(£)N(£) ~ V(£)N(£) 

K ) ln(l/7ei c ) - 1 + ftic j ln(l/ft lc )-l + ft lc ~ |lnft lc | 1 j 

(lMerrittli2013i . equations (6.59)-(6.62)), i.e. 

F EhC (£) « toWpFLC^' (56) 

| In 7vi c | 


We assume that equation (1511) describes F{£) for £ < £ CI - lt and that equation (1551) describes F{£) for £ > f C rit, 
where f C rit is the energy separating the full- and empty-loss-cone regimes. Identifying /(£), N(£) and T>(£) with 
the expressions for an unmodified Bahcall-Wolf cusp, as given in (J3l the total loss rates from the two regimes can be 
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written in terms of £„it after integration over £, as: 


N flc 

jyELC 


p 2 r(ll/ 4 ) M m ( GM.f 2 r lc f £ CI it \ 5/4 
V 7r r(5/4) TO* c 2 r ^/ 2 r g \ / 


25-\/2 

“32” 


C 


rr(n/4)i 

2 In A 

IGM . 

{ fcrit ^ 

- 1 c _ GM. 

[ r(5/4) J 

ln^V 

r 3 

1 m 

\£m) 

i C’m — 

^ m 


(57a) 

(57b) 


The constant C ~ 1.62 is defined in equation (1531) . In the integral for N ELC ; the term ln77i c was assumed independent 
of £; a reasonable choice for this term might be its value at £ = £ C rit, he. In77(£ cr it). 



m */ M 0 

Fig. 8.— Loss rates for the models of Figure[7] Scaling assumes M. = 4 x 10 6 M©. Plotted points are the same numbers given in Table fTI 
dashed and dotted curves are the approximate analytic relations derived in the text for the full- and empty-loss-cone regimes respectively, 
and their sum is shown as the solid curves. The open circles are from numerical integrations in which the £- grid probably did not extend 
to low enough values to yield the correct, total loss rates (see the right panel of Figure fTI). 


TABLE 1 

Steady-state loss rates 


P(r = 1 pc) 

m*/M 0 

N (jn 

f- 1 ) 

7V elu /A 

(M 0 pc- 3 ) 


(numerical) 

(analytic) 

1.9 x 10 4 

0.3 

(9.03 

X 

10 -b ) 

0.65 


1.0 

(7.34 

X 

10" 6 ) 

0.62 


3.0 

(5.74 

X 

10"®) 

0.60 


10. 

4.24 

X 

10“® 

0.58 


30. 

3.14 

X 

10“® 

0.55 


100. 

2.20 

X 

io-® 

0.53 

3.5 x 10 s 

0.3 

1.06 

X 

10“ 3 

0.59 


1.0 

7.12 

X 

10" 4 

0.56 


3.0 

4.91 

X 

10" 4 

0.54 


10. 

3.29 

X 

10" 4 

0.52 


30. 

2.31 

X 

10" 4 

0.50 


100. 

1.58 

X 

10 -4 

0.47 

6.1 x 10 6 

0.3 

1.03 

X 

10" 1 

0.53 


1.0 

6.45 

X 

10 -2 

0.51 


3.0 

4.22 

X 

10 -2 

0.48 


10. 

2.69 

X 

10 -2 

0.46 


30. 

1.81 

X 

10 -2 

0.44 


100. 

1.18 

X 

10 -2 

0.42 
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We take for £ cr j t the energy that satisfies 

9ic(£crit) = | lnft-id (58) 

since at this energy, equation (l56l) suggests that F ELC « F flc . Using equations (H|, (1201) and (l34l) we find 

_ 5 v ^r(ii/4) / 


«c(£) = ( 7 — 

v £"m 


£ \ °^ 4 to* In A 


M.n, c {£Y qo 2 r(5/4) 

When solving for £ cr it £\ c , the expression for q\ c can be simplified by writing 1Zi c (£) ~ 2£/£\ c , or 

,c x go to* £\ c (fcrit \~ 9/4 

Equating this with — lnTvhc yields a transcendental equation for x = £ CI it/£\c'- 


-<7« 12.73. 


cc 9 / 4 In (2x) = — — In A ( 


2 M. 


c \ 5/4 
Orr 




£\cj 


(59) 


(60) 


For values of x in the range of interest (10 7 < x < 10 5 ), an approximate solution to y = a; 9 / 4 ln(2x) is x = \J—2y, 
so that 


^crit 

2.75 

/to* 

£rn 

00 

0 

V M. 

lnP]- 1 * 

a-In 

11.0 0 


3/8 


M. 


5/8' 


(61a) 

(61b) 


Equations m and ED are the desired expressions. The predicted values for N FLC and /V ELC are plotted as the 
curves in Figure [8] The agreement is quite good considering the approximations made; one implication is that the 
modifications to the Bahcall-Wolf cusp resulting from resonant relaxation have little effect on the total loss rate. The 
predicted ratios N Fhc / (_/V ELC + fV ELC ) are given in the final column of Table [D in all cases considered here, the two 
regimes contribute roughly equally to the total loss rate. We emphasize again that these results apply only to the 
model considered here, which does not include the contribution of the distributed (stellar) mass to the gravitational 
potential. 

Another interesting property of the steady-state models is the distribution of orbital elements of the captured stars. 
In the empty-loss-cone regime (qi c <C 1), loss-cone orbits will have orbital periapsides close to n c , the physical radius 
of the loss sphere, but in the full-loss-cone regime (q\ c 1), timescales for change in L are comparable with orbit al 
perio ds and stars at the time of capture can be on orbits with every value of r p from zero to n c dCohn fe Kulsrudl 
1978). O ne consequence is that stars can e xperience stronger tidal stresses than if they were all lost from orbits with 
r p = 7i c (lOuillochon fc Ramirez-RuI3[2013 ). 

Figure [9] shows d 2 F/dr a dr p , the contribution to the loss-cone flux from stars with orbital apoapsides in the range 
r a to r a + dr a and orbital periapsides in the range r p to r p + dr p , for three steady-state models. The details of the 
calculation are given in the Appendix, which also presents the results of a similar calculation for the classical Bahcall- 
Wolf solution. As shown there, the computed distribution depends on two quantities: qi c {£), and /(£, 1Zi c ) = fic.(£)- In 
the Bahcall-Wolf model, unmodified by resonant relaxation, Figure fl2l shows that the distribution of orbital elements 
is a strong function of binding energy. At high £, i.e. small r a , q\ c 1 and the distribution of periapsides is strongly 
peaked toward r p = n c . At low £, i.e. large r a , the loss cone is full and F is nearly independent of r p 0 After 
integration with respect to r a , the distribution of orbital periapsides in the classical Bahcall-Wolf solution has the 
form shown in the right-hand panel of Figure fl2l with a maximum at r p = ?’i c and a very steep drop for r p < r\ c . 

As Figure [9] shows, the inclusion of resonant relaxation modifies this distribution, in the sense of reducing the 
contribution from stars with r p ss r\ c . This is an indirect consequence of the strong depletion in / at large binding 
energies. The depleted orbits are mostly in the enrpty-loss-cone regime, and their removal implies a larger relative 
contribution to the loss-cone flux from stars in the full-loss-cone regime, hence r p < r\ c . Figure [9] shows that this 
modification is severest in the model with the largest r m , i.e. the largest core; in this model, the distribution of 
captured stars with respect to r p is nearly uniform. Even in the model with smallest r m , i.e. the smallest core, the 
distribution with respect to r p is only mildly peaked near r\ Cl much less so than in the classical solution. 

The distributions shown in Figure [9] are computed from models with to* = M 0 . Models with larger to* have larger 
gi c , and the distribution of captured stars with respect to periapsis in these models is even more uniform than shown 
in Figure [9] 

To a reasonable approximation, therefore, one can write for all these models: 

P(<r p )«^, (62) 

He 


3 These properties of the orbital distribution were noted previously by IStrubbel (120111 ). 
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Fig. 9. — Distribution of orbital elements of stars fed to the SBH, computed as described in the A ppendix, for the three steady-state 
models from the right-hand panel of Figure [4] Left: Joint distribution of r a and r p , equation dAlOl ). Contours are spaced uniformly in 
log F and the range in contour values is 10 3 . Right: Thin curves show dN/dr p , obtained by integrating the function in the left panel with 
respect to r a at each r p . Thick curves are the number of stars with periapsides less than r p , obtained by a second integration with respect 
to r p . Both distributions are normalized to unit total number of stars. 

where P is the probability of capture from an orbit with periapsis less than r p . 
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4.2. Time-dependent solutions 

In their pioneering work, Cohn & Kulsrud (1978) presented only steady-state solutions. This was in keeping with 
their parameter choices, which were appropriate to massive black holes at the centers of globular clusters. But energy 
diffusion timescales near the centers of galaxies are often much longer than in globular clusters, and it is likely that 
many nuclei have not yet reached steady states under the influence of gravitational encounters. This is presumably the 
case in all galaxies with parsec-scale c ores; but even in galaxies with dense, nuclear star clusters, inferred relaxation 
times are often of order 10 9 yr or more (lMerrittll2009S) . The nucleus of the Milky Way probably falls in the non-relaxed 
category (IMerrittH2010f) . 

The time-evolution of such nuclei will differ depending on their assumed initial state. One widely discussed model 
invokes a binary SBH, which scatters and redistribu tes stars bef ore (presu mably) achieving a small enough separation 
that coalescence of the two black holes can occur dBegelman et al.lll980f). The late evolution of such binaries is not 
well understood, but their initial evolution appears to be fairly robust (TMilosavlievic fe Merrittl [20011) . After forming 
a bound pair, at a separation roughly equal to the influence radius of the larger SBH, the binary separation rapidly 
decreases through the combined influence of dynamical friction and three-body interactions with stars. This phase 
ends at the so-called “hard binary” separation, a ~ ah, where 


= G/i _ M 2 r h 
h 4 a 2 M\2 4 


(63) 


(lMerrittll2013L equation 8.23). Here, M i2 = Mi + M 2 is the binary mass; /r = M 1 M 2 /M 12 ; and = GMi/a 2 is the 
influence radius of the larger SBH. At separations < ah, the binary is able to eject stars with high enough velocities 
that they escape from the nucleus. The binary may “stall” at this radius; or, if a mechanism exists for repopulating 
the depleted orbits, its semimajor axis can continue to drop. 

During the early, rapid phase of its evolution, the binary interacts with stars on orbits having periapses (defined 
with respect to the binary center of mass) from ~ rh to ~ ah or less. This interaction modifies orbits with a range of 
periapses, from ~ down to ~ ah- Stars on orbits with initial periapses < ah are removed entirely from the nucleus. 
Here, we approximate the stellar distribution at the end of this phase simply as 


f{E,L) = f 0 (E), T>i g ap 

= 0; L -^gap; 


(64) 


with _ 

Lgap(E) = Ka h y/2[E-^{Ka h )], Knl, (65) 


the angular momentum of an orbit with periapsis at Kah■ The corresponding density profile has a core of radius ~ a 
few xah- 

Unless the binary mass ratio is close to unity, the core will be small compared with rh- In this circumstance, one 
expects post-binary evolution of the stellar distribution to take place on two timescales. Initially, the gap in / at low 
angular momenta is refilled via diffusion in L. The associated timescale is 

V-\£ gap ) (66) 


with f gap w GM,/ah- After a time of ~ T gap , the phase-space density in the 
will be approximately constant with respect to L at each E. 

On longer timescales, of order T gap < A t < T r . the distribution of orbital 
the Bahcall-Wolf steady state. 

Evolution of / in the first phase can be approximated by ignoring energy 
equation at each energy as 


dN 

~dt 




region previously emptied by the binary 
energies will evolve, eventually reaching 
diffusion and writing the Fokker-Planck 


(67) 


dMilosavlievic fc Merrittll2003 1. By changing variables from 1Z to £ = VE., equation (l67l) becomes the heat conduction 
equation in cylin dric al coor dinates with radial variable t and diffusivity T> and has a known solution in terms of 
basis functions dMerrittlI2013L 6.1.5). This solution has been applied to galactic nuclei, assuming for V the classical 
angular-momentum diffusion rate, i.e. 'D~ 1 ss T r (IMerritt fe Wang||2005f) . 

Initial conditions for the integrations presented here were constructed in the same way as in the previous section, 
with the added step of setting / to zero at L < L gap (U). Equation (1551) is not quite appropriate here given that ah is 
defined in terms of cr, the stellar velocity dispersion beyond the SBH influence sphere. Instead, the maximum periapsis 
of an evacuated orbit was computed from the roughly equivalent expression 


M r m 

fp.max - 4 


q r m 
(1 + q)' 2 4 


( 68 ) 


with q = M 2 /M\ the binary mass ratio (IMerritt fe SzeHlf2006! ). 
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TABLE 2 

Evacuated-core models 


7 

<? 

^p,max (pc) 

3/2 

0.01 

0.005 


0.03 

0.014 


0.1 

0.041 

2 

0.01 

0.015 


0.03 

0.042 


0.1 

0.124 


Six initial models were constructed: one set with 7 = 3/2 and one set with 7 = 2. These values are, respectively, 
less than and greater than the Bahcall-Wolf value 7 = 7/4. For each 7, the binary mass ratio in equation (1551) was 
assigned one of the three values 

q = {0.01,0.03,0.1} (69) 

and the initial / was set to zero for orbits with periapses less than the r Pjmax given by equation (1681) . with K = 1. The 
stellar mass was set to 1M 0 , assuming a SBH mass of 4 x 1O 6 M 0 , and r\ c was set to 15 r g as in the previous section. 
The other important initial parameter was the value of the mass density at large radii. This was set to 

p(r = 1 pc) ~ 4.0 x lO 5 M 0 pc -3 

in the three models with 7 = 3/2, and 

p(r = 1 pc) ~ 1.0 x lO 5 M 0 pc -3 

in the models with 7 = 2. The corresponding values of r m were approximately 2.0 pc and 6.0 pc, respectively. Tablc[2] 
lists the important parameters for the six initial models. 




Fig. 10.— Evolution of p(r) in two integrations starting from evacuated-core initial conditions, with q = 0.03. Left panel: 7 = 3/2; 
right panel: 7 = 2 . The initial density profile is shown as the dashed curve; subsequent times are displayed as curves of increasing width. 
Displayed times are {0.003, 0.01, 0.03, 0.1, 0.5,10} x 10 9 yr (left panel) and {0.01, 0.03, 0.1, 0.2, 0.5, 10} X 10 9 yr (right panel). Scaling assumes 
M m = 4 X 10 6 Mq. In the right-hand panel, note that the final central density is slightly lower than its value at t = 5 X 10 s yr. 


Figure[10]shows the evolution of the mass density profile in the two integrations with q = 0.03. The values of r p ^ mSLX 
were ~ 0.014 pc (7 = 3/2) and ~ 0.042 pc (7 = 2). The two evolutionary timescales discussed above are evident. 
After a time of ~ 10 8 yr (7 = 3/2) or ~ 3 x 10 8 yr (7 = 2), the initial core is erased; during this time, the density at 
r v p max hardly changes. Over the next 10 9 — 10 10 yr, energy relaxation causes p(r) to approach the Bahcall-Wolf 
form at large radii. 
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Evolutionary models like these are characterized by an additional dimensionless parameter: the ratio of the initial 
core size, ~ r Pimax , to the size of the core that forms via resonant relaxation, in the manner discussed in the previous 
section. Adopting the estimate given above, r c « 0.03r m , for the latter core size, the ratio becomes 


•p, max q Q 


(70) 


This ratio is ~ 0.25 for both of the models of Figure [TUI implying that the final core should be somewhat larger than 
the initial core, as seen in the figure. There is however a change in the structure of the core. Initially, the core is 
formed by the exclusion of orbits with small periapses, implying a zero configuration-space density below some radius. 
The core that forms at late times is characterized by a deficit of orbits with high binding energies; as discussed above, 
the implied density is nonzero, p ~ r -1 / 2 near the center, due to orbits with low binding energies and small angular 
momenta that pass near the center. 

The fact that the ratio in equation (ITUl) is less than unity for these models implies that the initial timescale for core 
refilling is set by resonant, and not classical, relaxation. We can estimate that time from equation (1361) . replacing a 
by r pmax . The result is ~ 6 x 10 7 yr (7 = 3/2) and ~3x 10 8 yr (7 = 2), quite consistent with the evolution seen in 
Figure EB 



time (yr) time (yr) 

Fig. 11.— Total loss rates for the six models with evacuated cores at t = 0. Left: 7 = 3/2; right: 7 = 2 . Scaling assumes M, = 4x 10 6 Mq. 
Dashed curves are from integrations starting from the same initial conditions, but with the energy diffusion terms artificially set to zero. 


The two evolutionary timescales are reflected also in the change with time of the SBH feeding rate N. Figure |TT) 
plots this quantity for each of the six models with initially evacuated cores. Also plotted there, as the dashed curves, 
are the feeding rates from a second set of integrations in which the energy diffusion terms were artificially set to zero. 
In those models, evolution of / at each E is described approximately by equation (El). Comparison of the two sets of 
curves shows in a very direct way how the early evolution is determined by angular-momentum diffusion and the late 
evolution by energy diffusion. 

In the models of Figure [Tl] steady-state loss rates are only achieved after a time of ~ 10° yr (7 = 3/2) or ~ 3 x 10 9 
yr (7 = 2). These times are fixed by the classical relaxation time, equation (j4j, and hence by the adopted density 
normalization. It is clear from this small set of examples that many nuclei, including that of the Milky Way, might not 
be in a steady state with regard to SBH feeding rates. Indeed, in nuclei with SBH masses larger than the value assumed 
here (4 x 10 6 Mq), equation (l36l) implies that even the initial timescale for core refilling due to angular momentum 
diffusion could exceed 10 10 yr. 


5. DISCUSSION 


5.1. Comparisons with earlier work 

The consequences of resonant relaxation for the steady-state distribut ion of stars around a SBH have been discussed 
by earlier authors (Hopman & Alexander 2006b; Madiga n et al.ll201lD using more approximate methods. Here we 
com pare the resul ts of th ose stu dies with the results obtained here. 

iHomnan fc Alexander! ( 2006bl i derived steady-state solutions for N ( E ) in the fixed gravitational potential of a SBH. 
Their evolution equation had the form 


d N 
~dt 


dF e 

dS 


■FnrOM) - X 


N 

Trr 


(71) 
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The ter m Fnr accounts for loss of stars into the SBH vi a “non-resonant,” i.e. c lassical, diffusion in Z/: lHopman fc A lexander 
lj2006bT ) adopted an expression similar to that used bv lBahcall fc WoH (jl977f ) in their earlier study of the Id problem. 
The last term on the right hand side of equation m approximates the loss rate at energy £ due to resonant relaxation; 
Xrr is an estimate of the resonant relaxation time, and the dimensionless factor \ — 0(1) was included to parametrize 
uncertainties in the efficiency of resonant relaxation and in the degree of depletion of phase space near the loss-cone 
boundary; the latter could not be modeled in their study due to the Id, / = f(£) approximation. 

We can cast our Fokker-Planck equation for /(£,77,7) in an analogous form, as follows. Equation (fTUl) is 


dN(E, 77) 

at 


<9 


<977 


df 


= ■ ■. J 7— Dim—— + Duf 


<977 


(72) 


where terms depending on the energy-space flux have been omitted. Adopting the expressions (1131) for the flux 
coefficients due to RR, equation (l72l) becomes 


dN (g.K) _ 2A(£) — 

dt {t ’dTZ 




Integrating this expression <777 yields an evolutionary equation for N(£): 


dN(£, t) 
dt 


= ... 2 A{£) 




- Wlc 


= • ■ • — 2A(£)77i c (£) [1 — 77i c (£)] 


dN \ 

an) 


(73) 

(74a) 

(74b) 


^lc 


To make further progress we need an ansatz for the 77-dependence of N. A natural choice is the “empty-loss-cone” 
solution (1531) . which implies 


77(77; £) = N(l;£) 1- 


ln77 
In 77i c 


N{£) In (77/77i c ) 


In 77, 


-i 


77i c — 1 


where we have identified N(£) = f 77(77; £}d77. Inserting (1751) into (I74bl) then yields 


dN 

dt 


and comparing this expression with dm we conclude 


2A(£)N(£) 

ln(l/77i c ) 


T RR (£)A(£) 

ln(l/77i c ) 

Adopting equation (fl4l) for A(£), and iHopman fc Alexander! (l2006bf) ’s expression for Trr: 


(75) 


(76) 


(77) 


= Arr fM.\ 2 P 2 (q) 
7V*(< a) V TO */ Aoh(a) ’ 


Arr ~ 3.56 


(78) 


then yields 
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l ^ 1 ' 


(79) 


In the models computed here, equation (I61bl) gives ln77 1 7 1 « 12, so that \ ~ 1.5. iHopman fc Alexander! (I2006bll 
presented steady-state solutions for the cases \ = {1, 3,10}. Their Figure 8 shows steady-state density profiles, p(r), 
for the Milky Way nucleus, assuming x = 1 and x = 10; the re are depletions wit h r espect to the x = 0 (classical) 
case inside radii of ~ 0.05 pc and ~ 0.2 pc, respectively. While IHopman fc Alexander! (l2006bT) do not clearly state the 
density normalization for their models (there are no units on the vertical axis of their Figure 8), the mass density at 
1 pc appears to be ~ 1 x 1O 5 M 0 pc -3 , which would place it midway between the two models in Fi gure |H w ith t he 
l ower d ensity normalizations. The core sizes in those models are similar to the values found bv IHopman fc Alexander! 
(I20061J) . 

One issue that complicates a comparison of Hopman & Alexander’s (2006) results with ours is the choice that they 
made for the coherence time. They defined 


1 

Aoi, 


1 


1 


^coh.M 


(80) 


with fcoh.M and f C oh,s defined in essentially the same way as here (equation fl5l) . Both quantities are positive by 
definition, and the minus sign in equation (1801) was said to account for the fact that mass precession is retrograde 
and Schwarzschild precession prograde. At a certain energy / radius, the coherence time as defined by equation (1801) 
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becomes infinite. While the precession rate of a single orbit can be zero (if its eccentricity has precisely the right 
value), orbits of other stars at similar radii will still precess, implying a finite_coherence time at every radius. 

A very different approach to the problem was taken bv iMadigan et all (120111) . who developed an ad-hoc statistical 
model for the effects of resonant relaxation, calibrated against IV-body simulations. Nevertheless, their and our 
results about the depletion of f(E) at high binding energies and the correspo nding flattening of the density profile 
are qualitatively similar. The size of the steady-state core predicted bv IMadigan et ah! (1 2011 ) for the Milky Way was 
< 0.05 pc, consistent with the core sizes in some of the models shown here in Figure © 

5.2. Comparison with the distribution of stars at the center of the Milky Way 

The distribution of stars near the center of the Milky Way has long been known to depart from the Bahcall-Wolf 
steady-state form. Solutions of the isotropic Fokker-Planck equation that include the stellar potential suggest that 
the Bahcall-Wolf solution (which ignores the stellar potential) should be valid out to distances of at least ~ 0.2 r m 
from the SBH ( IMerrittll2013L 7.1.1). Applied to the Milky Way, this result implies that the n ~ r -7 / 4 cusp, if present, 
would extend outward to ~ 0.2 — 0.6 pc. However, number counts of the late-type (i.e. old) stars fail to show a cusp. 
Instead, the density of these stars (most of whic h are believed to be red giants) rises only very slowly, if a t all, toward 
t he center inside a proj ected radius of ~ 0.5 pc (|Buchholz et alJl2009HDo et al.ll2009l : iBartko et al.ll2010f) . 

IMadigan et all (1201 lh proposed that resonant relaxation was responsible for the Milky Way core. However the 
steady-state cores found by th ose auth ors were about an ord er of magnitude smaller than the core observed in the 
Milky Way. The core sizes in iHonman fe Alexander! (|2006bD ’s steady-state models were also substantially smaller 
than 0.5 pc, as noted above. The arguments presented here (§ © suggest that the size of a steady-state core that is 
produced by the action of resonant relaxation should scale roughly with r m . Based on the results plotted in Figure [2 
a steady-state core as large as ~ 0.5 pc in the Milky Way would require r m > 20 pc, far larger than most estimates of 
the SBH influence radius. 

A likely resolution to this apparent di screp anc y is t o assume that the Milky Way nucleus is not yet in a steady 
state with regard to diffusion in energy (|Merrittl(2010f ). The observed core could then result from a combination of 
initial conditions, which have not yet been erased by the effects of gravitational encou nters; and t he depleting effects 
of resonant relaxation. This possibility is explored in more detail in a separate paper (|Merritt et alJl2015h . 

5.3. Other consequences of the depletion of f at large binding energies 

The striking depletion in N(E) at large binding energies found here, and in some earlier studies, is a consequence 
of two things: (i) the sudden drop in the angular momentum diffusion time below a certain energy, due to resonant 
relaxation; and (ii) the absence of any mechanism that might maintain a high density of stars near the SBH in spite of 
the high loss rate. The first assumption seems robust; the second less so, since the forms for the diffusion coeffcients 
adopted here are only likely to be valid beyond a certain distance from the SBH. If diffusion times become long again in 
the “Schwarzschild,” “Kerr” etc. regions of Figure 1, a high density of stars (or compact objects) might be maintained 
very near the SBH. Ongoing star formation in this region could also maintain a nonzero /. 

Assuming for the moment that neither happens, one can speculate about the consequences of the depletion. An 
approximate representation of the steady-state phase space density in the models of Figures ©and [5] is 

/(£, ft) «/(£), £<£ eq 

wo, £>£ eq (81) 

with £ eq the binding energy at which resonant relaxation begins to dominate classical relaxation roughly, the energy 
of an orbit wit h radius a fin given by equation (1551). Pr operties of models having the functional form (I81|) for / have 
been discussed ((Me rritt 20l0[ lAntonini fe Merrittll2012f) . The distribution of orbital eccentricities, N(e), approximates 
a delta function, N(e) ~ 5(1 — e), for r <C a eq , since the only orbits that approach closely to the SBH are very eccentric 
ones. 

Figure © suggests that in the Milky Way, the region of strong depletion would have a radius < 0.1 pc. As noted 
above, this is probably too small to explain the observed core; but since the depletion due to resonant relaxation 
occurs ona~ 10 8 yr timescale, it is likely to be present whatever the explanation for the larger observed core unless 
the “initial conditions” were extreme. 

In the Milky Way, the brightest stars in this region are the S-stars and the stars in the two stellar disks. Their 
presence is not inconsistent with the depletion discussed here since these stars must have formed very recently: less 
than ~ 10 8 years ago in the case of the S-stars, and less than ~ 10 7 years ago in the case of the stellar disks dSchodell 

unnD. 

The evolutionary models presented here are more relevant to old stellar populations. Stellar-mass black holes, 
with masses ~ 3 Mp — 30are likely to domin ate the mass density inside ~ 10 ~ 2 pc from the Milky Way SBH 
((Freitag et all 2006: H ooman fc Alexander! I2006al) . due to mass segregation and due to ti dal destruction of norma l 
stars. These are the objects that could become EMRIs, or extreme-mass-ratio inspirals ((Sigurdsson fc Re© 1 19971) . 
of great interest to experimental physicists hoping to detect low-frequency gravitational waves. The depletion in 
/ discussed here is likely to have important consequences for the steady-state rate of EMRI production. Exactly 
what those consequences are can not be stated with certainty yet, since the very eccentric orbits that lead to EMRIs 
probably evolve in a way that is not well described by the low-L forms of the diffusion coefficients assumed here. 
Rather, these orbits are subject to “anomalous relaxation,” the qualitatively different way in which orbits evolve when 
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their precession rate (due to GR) is much higher than that of the field stars (|Merritt et al.l[20TTI : lAntonini fe Merritt! 

EoH. 

6. SUMMARY 

Integrations of the Fokker-Planck equation describing /(E, L,f), the phase-space density of stars around a super- 
massive b lack h ole £SBH) at the center of a galaxy, were carried out using a numerical algorithm described in an earlier 
paper (iMerrittll2015f ). Diffusion coefficients describing both classical and “resonant” relaxation were included. Both 
steady-state and time-dependent solutions were found. The principal results follow. 

1. Steady-state solutions, with fixed density far from the SBH, are similar to the classical, isotropic, Bahcall-Wolf 
solution, i.e. / ~ \E \ 1 / 4 , n ~ r 7 / 4 . However the enhancement of angular momentum diffusion at large binding 
energies, due to resonant relaxation, implies a depletion in / at those energies and a corresponding density deficit, or 
“core.” The core radius scales approximately with the gravitational influence radius of the SBH and is a few percent 
of that radius. The density within the core is n ~ r~ 1 ^ 2 . 

2. Although the inclusion of resonant relaxation has a substantial effect on the density profile near the SBH, the 
consequences for the SBH feeding rate are much less extreme, since most stars are scattered into the SBH from orbits 
that lie outside the resonant relaxation regime and since / is strongly depleted in that region. A simple analytic 
formula is derived, based on the classical diffusion coefficients, that reproduces the numerically-computed loss rates 
with good accuracy. The depletion in / at large binding energies does significantly affect the distribution of orbital 
elements of captured stars, in the sense of reducing the contribution from stars on orbits with periapsides r p ss n c , the 
radius of the physical loss sphere. 

3. Since energy relaxation times at the centers of galaxies are often very long, time-dependent solutions were also 
computed. Initial conditions were based on a model in which the SBH was preceded by a massive binary. These 
models evolve on two timescales: a short timescale during which the core evacuated by the massive binary is refilled 
via angular-momentum diffusion; and a longer timescale during which diffusion in energy causes the radial distribution 
of stars to approach the Bahcall-Wolf form far from the SBH. 

4. Steady-state cores produced by the effects of resonant relaxation in these models are probably too small to explain 
the core observed in the distribution of late-type stars at the center of the Milky Way. A possible resolution of the 
apparent discrepancy would be to assume that the nucleus of the Milky Way has not yet reached a steady state under 
the influence of gravitational encounters. 

5. The depletion in / at large binding energies could have important consequences for the production of EMRIs, or 
extreme-mass-ratio inspirals. A final decision on this question must await a more careful treatment that includes the 
effects of “anomalous relaxation,” the qualitatively different way in which eccentric orbits evolve in the post-Newtonian 
regime. 


This work was supported by the National Science Foundation under grant no. AST 1211602 and by the National 
Aeronautics and Space Administration under grant no. NNX13AG92G. 

APPENDIX 

DISTRIBUTION OF PERIAPSIDES 

In the Cohn-Kulsrud boundary layer treatment, at the moment of capture by the SBH, stars can be on orbits with 
periapsides in the range 0 < r p < HcD The distribution of orbital integrals {£, TV) of captured stars is given by equation 
(6.57) of lMerrittl (|2013f ) (hereinafter DEGN). Here we recast that equation in terms of orbital elements, including the 
periapsis distance r p , and use the result to compute the periapsis distribution for stars in a classical Bahcall-Wolf cusp. 
The same algorithm was used to compute the periapsis distribution for the Fokker-Planck models in § 14.11 iFigurelfll). 

Differentiating equation (6.57) of DEGN with respect to TZ yields the contribution to the loss cone flux from stars 
in the energy interval £ to £ + d£ and angular momentum interval TZ to TZ + cTJZ: 


d 2 F 

d£dTZ 


= 4n 2 L 2 (£)f(£,TL,T = l), TZ<lZ lc {£). 


(Al) 


In equation (IA1I) . the phas e space dens ity / is understood to be a function of position along an orbit, for any orbit 
having periapsis inside r\ c (iGohn fc Kulsrudl 119781) . The dimensionless variable r measures position along such an 
orbit; as r varies from 0 to 1, r increases from periapsis (r = r p , f = 0), to apoapsis (r = r a ) and back to r p again 
(where / has its maximum value along the orbit). The quantity f(£,TZ,r = 1) is expressible in series form as 


f(£,TZ,l) = f(£,TZ lc ) 


1 - 


2 y e~P 2 m/* Jo(PmVy) 

\/2/lc m=1 Pm Jl{Pmy/V\c 


(DEGN, equation 6.54). In equation 


y is a dimensionless angular momentum variable defined as 
TZ TZ\ C 


V = 


P(£)V(£) ’ Vlc P(£)V{£) ’ 

4 The calculation described here is non-relativistic and no special consideration is given to stars that travel near or 
inside the hole’s event horizon, r < r g . 
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(A3) 
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Jo and J\ are Bessel functions of the first kind, and the /3 m yield successive zeros of the equation 

Jo {Py/Vte) = 0 . 

In terms of a m = ftmy/yic and x = y/y q\ c , equation (IA2I) is 

f(£,U,l) = f lc (£)W(q lc ,x), 

jxn \ 1 0 W' e _Q ”* 91c / 4 J 0 (a m x) 

W(qic, x) = l-2\ - — - r 

, a m J i (a m ) 

m= 1 v 7 

where qi c = q\ c (£) = t/jy (£) is defined as in equation (EU1) and f\ c (£) = f[£, 1Z\ C (£)\. 

In the Kepler potential assumed here, L 2 (£) = G 2 M, 2 / (2£) so that 

= 2n 2 (GM.) 2 £- 1 f lc (£)W(q lc , x). 

This can be converted into a distribution in (r a , r p ) using 

r a = a (1 + e) = (l + Vl - Tlj , r p = a (1 - e) = (l ~ Vl -II s ) 


i.e. 


The Jacobian is 


£ = 


GM. 


n = 


4 r a r p 


d(£,H) _ 4£ 3 


( r a + r p ) 


;VTJJtz = 4GM.- 


2 • 


'a 1 p 


so that 


d(r a ,r p ) G 2 M, 2 
d2p -, 2 ^ 21*2 (r a -r p ) 


dr a dr p 


= 8 ? t 2 G 2 M. 


(■ r a + r p ) 3 


(r a + r p ) 

,f\c{£)W(q\ c ,x). 


(A4) 

(A5a) 

(A5b) 


(A 6 ) 

(A7) 

(AS) 

(A9) 

(A10) 


The quantity x = x(£, TV) = y/TZ/TZi c {£) that appears in these expressions can be written in terms of the orbital 
elements as 


x = 


' a' p 


n C (r a +r p -n c ) ' 


(AH) 


In the limit of nearly-unbound orbits, i.e. £ ~ 0, x ~ y/r p /r\ c . 

These relations can be applied to a classical Balicall-Wolf cusp. The function q \ c (£) follows from equations (EU1) . 
( 1 M 1 ) and (fTTTTl] : 

q\c{£) = 8/ 0 \/27r 3 lnA7e i - 1 (f)(GM.) V4 Gm*r- 5 / 4 £- 5 / 4 , f 0 « 0.036. (A12) 


We relate fic(£) to f(£) using equations (6.61) and (6.63) of DEGN: 


U£) 


m 

1 + « , lc 1 ^(®c) ln ( 1 /^lc) 


(A13) 


and identify / with equation (l26l) . The results are shown in Figure HU for a nucleus with r m = 10 7 r g , m* = 
M./(4 x 10 6 Mq), InA = 15, and r\ c /r g = 15. The computation included 5000 terms in the Bessel series and used a 
(500 x 500) grid in (r a , r p ). The distribution with respect to r p can be seen to become nearly uniform for large r a (full 
loss cone), while for small r a (empty loss cone), the distribution is very strongly peaked toward r p = ri c . Roughly 1/2 
of all captured stars have periapsides in the range 10r a <r p < 15r ff (= n c ) (right panel). 
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Fig. 12.— Distribution of orbital elements of stars fed to the SBH in a classical Bahcall-Wolf cusp, according to the Cohn-Kulsrud 
boundary-layer solution. Parameters were r m = 10 7 r^, m* = M*/(4 x 10 Q Mq), In A = 15, r\ c /r g = 15. Left: Joint distribution of r a and 
r p , equation (IaToIi . Contours are spaced uniformly in logF and the range in contour values is 10 3 . Right: Thin curve shows dN/dr p , 
obtained by integrating the function in the left panel with respect to r a at each r p . Thick curve is the number of stars with periapsides 
less than r p , obtained by a second integration with respect to r p . Both curves are normalized assuming a unit total number of stars. 
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